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Magnetic Reconnection in a Weakly Ionized Plasma 
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Magnetic reconnection in partially ionized plasmas is a ubiquitous phenomenon spanning the range from lab- 
oratory to intergalactic scales, yet it remains poorly understood and relatively little studied. Here, we present 
results from a self-consistent multi-fluid simulation of magnetic reconnection in a weakly ionized reacting 
plasma with a particular focus on the parameter regime of the solar chromosphere. The numerical model in- 
cludes collisional transport, interaction and reactions between the species, and optically thin radiative losses. 
This model improves upon our previous work in Leake et al. 2012^^ by considering realistic chromospheric 
transport coefficients, and by solving a generalized Ohm's law that accounts for finite ion-inertia and electron- 
neutral drag. We find that during the two dimensional reconnection of a Harris current sheet with an initial 
width larger than the neutral-ion collisional coupling scale, the current sheet thins until its width becomes less 
than this coupling scale, and the neutral and ion fiuids decouple upstream from the reconnection site. During 
this process of decouphng, we observe reconnection faster than the single-fluid Sweet-Parker prediction, with 
recombination and plasma outflow both playing a role in determining the reconnection rate. As the current 
sheet thins further and elongates it becomes unstable to the secondary tearing instability, and plasmoids are 
seen. The reconnection rate, outflows and plasmoids observed in this simulation provide evidence that mag- 
netic reconnection in the chromosphere could be responsible for jet-like transient phenomena such as spicules 
and chromospheric jets. 



I. INTRODUCTION 

Magnetic reconnection, the process which converts en- 
ergy stored in the magnetic fleld into kinetic and thermal 
energy of the plasma by breaking magnetic connectiv- 
ity, is observed in a range of laboratory and astrophysi- 
cal plasmas. In the astrophysical context, these include 
the interstellar medium (ISM), galactic disks, and the 
heliosphere^. In the heliosphere in particular, magnetic 
reconnection is believed to be the cause of transient phe- 
nomena such as solar flares and X-ray jets^, magneto- 
spheric substorms^^^ and solar 7-ray bursts^^^. 

Reconnection occurs when the ideal MHD frozen-in 
constraint, valid for highly conducting plasmas, is bro- 
ken locally, by allowing fleldlines to reconnect through a 
narrow diffusion region. Parker and Sweet (Refs. and 
[1) were the flrst to formulate magnetic reconnection as a 
local process within the single-fluid, fully ionized, magne- 
tohydrodynamic (MHD) framework. They considered a 
current layer of width Ssp much smaller than its length L, 
and plasma with non-zero resistivity {rj) due to electron- 
ion collisions. The model assumes steady state, i.e., the 
rate of supply of ions to the reconnection region by in- 
flow is equal to the rate of removal of ions by outflow, 
and takes the length L to be the characteristic system 
length scale. Using Ohm's law for a fully ionized, single 
fluid plasma with velocity v, magnetic fleld B, electric 
fleld E, and current density j: 



along with the plasma momentum equation and the 
steady state continuity equation, a simple equation for 
the reconnection rate can be derived: 
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Here S = iiqvaL/t] is the Lundquist number, Vin is 
the ion inflow upstream from the current sheet, va = 
B/y/pjjLQ is a typical Alfven velocity, with B evaluated 
upstream from the current sheet and p (mass density) 
evaluated in the current sheet. The permeability of free 
space is denoted by jiQ. From this analysis, it can also 
be shown that the Sweet-Parker current width 5sp and 
aspect ratio crsp = Ssp/L approximate to 
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E + vxB = 7^j, 



(1) 
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In this paper, we will focus on magnetic reconnec- 
tion in the solar chromosphere, which exhibits localized, 
transient outflows on a number of different length-scales. 
"Chromospheric jets" are the largest class of these out- 
flows, and are surges of plasma with typical lifetimes of 
200-1000 s, lengths of 5 Mm, and velocities at their base 
of 10 km/s (e.g., Refs. Will). "Spicules" are a smaller 
class of chromospheric outflows and have lifetimes of 10- 
600 s, lengths of up to 1 Mm and velocities of 20-150 
km/s (e.g., Ref. fll). A unifled model for the generation 
of chromospheric outflows was suggested by Refs. H and 
[9I. In their model, a bipolar fleld emerging into and recon- 
necting with a unipolar fleld creates outwardly directed 
reconnection jets. Depending on the size of the bipole 
and the strength of the fleld, this model may be able to 
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explain some of the observed properties of spicules and 
chromospheric jets. In Leake et al. 2012^^ we argued 
that, given the observed lifetimes, lengths and velocities 
of chromospheric jets, the minimum normalized recon- 
nection rate is M ^ 0.5. For spicules the minimum value 
we derived is M ^ 0.01. 

The solar chromosphere is a weakly ionized plasma, 
where the average collision time between neutrals and 
ions is on the order of ms, and so these two fluids can 
generally be considered to act as a single fluid. In this 
regime, the Sweet-Parker reconnection scaling is applica- 
ble, though the Alfven speed, va^ now depends on the 
total (ion+neutral) mass. However, if the width of the 
current sheet is comparable to the neutral-ion collisional 
mean free path, A^^ = VT^n/'^ni (where vt^u is the ther- 
mal velocity of neutrals and Vni is the neutral-ion col- 
lisional frequency), ions and neutrals can decouple and 
the mix may not be in ionization balance. In this case 
the reconnection region will have sources and sinks of 
ions in the form of the atomic reactions of ionization and 
recombination, and the inflow and outflow rates for the 
reconnection region could be affected by these interac- 
tions. 

The model for the multi-fluid simulations we presented 
in Leake et al. 2012^^ included ionization, recombination 
and scattering collisions, but did not include charge ex- 
change collisions. We found fast reconnection occurred 
when the ions and neutrals became decoupled upstream 
from the reconnection site, and the resulting recombina- 
tion of excess ions within the current sheet was as efficient 
at removing ions as the expulsion by reconnecting field. 
In this paper we improve upon that model in a number 
of ways. Firstly, we include charge exchange collisions 
and use more realistic cross-sections for scattering colli- 
sions, which effectively increases the collisional coupling 
of neutrals and ions, and thus decreases the neutral-ion 
collisional mean free path. At the same time, we also have 
a higher Lundquist number (see below) than in Leake et 
al. 2012i^ which reduces the Sweet-Parker width. As 
we will show, these changes result in similar ratios of the 
Sweet-Parker width to the neutral-ion collisional scale. 
We suggest that it is this ratio that is important for ob- 
taining fast reconnection, and so the simulation in this 
paper and those in Leake et al. 2012^^ are consistent 
despite these changes. Secondly, we allow the neutral 
fluid collisional transport coefficients to be functions of 
the local plasma conditions, rather than model parame- 
ters, and use more realistic values for the fluid transport 
coefficients. This results in a magnetic Prandtl number, 
the ratio of viscous to resistive diffusivity, which is lower 
than the simulations of Leake et al. 2012^^. 

Thirdly, we solve a generalized Ohm's law that ac- 
counts for finite ion inertia and electron-neutral drag ef- 
fects. In the simulations of Leake et al. 2012^^ the resis- 
tivity was a parameter of the simulations, and we derived 
scaling relationships for the reconnection rate and cur- 
rent sheet width in terms of the resistivity (or Lundquist 
number). The resistivity used in this paper, which ap- 



propriately depends on the plasma density and temper- 
ature, is a more realistic estimate for the chromospheric 
resistivity, and hence Lundquist number. The numerical 
model used in this paper is presented in §2 and the re- 
sults are presented in §3. §4 provides discussion of the 
results in reference to the transient phenomena observed 
in the chromosphere such as spicules and jets. 

II. NUMERICAL METHOD 

A. Multi-Fluid Partially Ionized Plasma Model 

The chromospheric model is based on the model used 
in our previous simulations in Leake et al. 2012^^, and 
so we only present here the modifications made to that 
model, and the necessary information to understand the 
modifications. The model consists of three fluids, ion 
(i), electron (e), and neutral (n). Each fluid {a) has a 
number density n^, particle mass pressure tensor 
Pa, velocity Vq,, and temperature Tq,. Only hydrogen 
is considered here, and so and are equal to the 
mass of a proton. These fluids can undergo recombina- 
tion, ionization and charge exchange interactions. The 
recombination and ionization reaction rates are given by 
Equations (4)-(7) in Leake et al. 2012^^. As mentioned, 
in Leake et al. 20 12^^ we did not include charge exchange 
interactions, but do so in the simulation presented in this 
paper. 

The charge exchange reaction rate F^^ is defined as 
r^- = Se,(K.)nin„V;,, (4) 

where 

Vex = \I^V^n + ^VTn+vl (5) 

is the representative speed of the interaction and vf^ = 
|vi — Vn|^. The thermal speed of species a is given by 
VTa = y^'^kBTa/rria, where /c^ is Boltzmann's constant. 
Following Ref. [l^, the charge exchange cross-section 
(Vex) is given by 

^cx{Vcx) = 1.12 X 10"^^ - 7.15 X IQ-^^lnVca: (6) 

which is a practical fit to observational data of Hydrogen 

collisionsi^. 

Continuity: 

The continuity equations for the ions and neutrals are 
given by Equations (10) and (11) in Leake et al. 2012—, 
respectively, and include ionization and radiative recom- 
bination. The electron continuity equation is not re- 
quired as charge neutrality (ng = Ui) is assumed. 
Momentum: 

The ionized fluid (electron+ion), and neutral fluid mo- 
mentum equations are given by Equation (12) and (13) 
in Leake et al. 2012^^, respectively, and include ioniza- 
tion, recombination and charge exchange contributions. 
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As in Leake et al. 2012— , we neglect the viscous part of 
the electron pressure tensor, but in this paper the con- 
tributions due to charge exchange are retained. 

In addition, we use different collisional cross-sections 
for ions and neutrals, compared to the model in Leake et 
al. 20 12^^. The collision frequency for collisions of ions 
with neutrals is is given by 



TTniin 



(7) 



. The ion neu- 



with T,^ = 2^ and m,;3 - • 
tral collisional cross-section is T^in- In this paper we use 
Tjin = Tjni = 1.16 X 10~^^ m^, whereas in Leake et al. 
2012^^, the model used = ^ni = 1-4 x 10"^^ m^ 
These improvements are made to take into account the 
theoretical work of Ref. 0. 

In the model in this paper we also use plasma de- 
pendent viscosity coefficients in the pressure tensor, un- 
like in Leake et al. 2012^^ where they were merely pa- 
rameters (the neutral viscosity coefficient was set to 
10~^ kg/(m.s), as was the ion viscosity coefficient). In 
this paper the neutral and ion viscosity coefficients are 
given by 

nnksTn - , ni^oksTio 
U = , and , (8) 

^nn ^ii 

respectively. The relevant collision frequencies Vnn and 
Uii are given by 



, and Uii = -niT^iix (9) 

Trmr, 3 V 7rm,- 



Ayrr^ ^, with 



where Snn = 7.73 x 10"^^ m^ and E^, 

= e^/(47reo/cBT^) being the distance of closest ap- 
proach for ions (eo is the permittivity of free space, e is 
the elementary charge, and A = 10 is the Coulomb loga- 
rithm). The subscripts in the ion coefficient represent 
background values at time t = 0. Hence the neutral vis- 
cosity coefficient varies but the ion viscosity coefficient 
has a constant value equivalent to unmagnetized ion vis- 
cosity in a uniform background plasma of temperature 
Ti^Q and density n^^o- Future work will consider the ef- 
fects of an anisotropic ion stress tensor. 
Internal Energy: 

The internal energy equations for the ionized fluid 
(ion-helectron) and the neutral fluid are given by Equa- 
tions (19) and (20) in Leake et al. 2012^^, respectively, 
and include ionization, recombination, charge exchange, 
and optically thin radiative losses. In this paper we re- 
tain the charge-exchange interaction contributions. We 
also use a neutral thermal conductivity {nn) that is de- 
pendent on the plasma parameters 



(10) 



This has a value of Kn = 1.76 x 10^^ m ^s ^ for the initial 
background neutral density and temperature used in the 



simulation in this paper. In Leake et al. 2012^ we used 
a constant value, independent of the neutral fluid condi- 
tions, of tZn = 1-58 X 10^^ m~^s~^. Note that the fric- 
tional heating and thermal transfer term was quoted 
incorrectly in Leake et al. 2012^^ and is instead given by 
Q-^ = ^R-/^ • (v^ - v«) + 3^naiya^kB{T^ - T«). 
Ohm's Law: 

The generalized Ohm's law used in this paper which the 
viscous part of the electron pressure tensor but includes 
electron-neutral collisions as well as electron-ion colli- 



sions: 



E + (v^ X B) = T^j ■ 



j xB 



-VPe 



-w, (11) 



where w = — v^. The resistivity rj is calculated using 
the electron-ion and electron-neutral collisional frequen- 
cies, and thus depends on the plasma conditions: 



V 



(ene)^ 



(12) 



where the electron- ion {uei) and electron- neutral (i^e 
collision frequencies are given by 



2kBTe 



and Up 



SksTe, 



(13) 

The electron- ion collisional cross-section Sgi = Ayrr^g, 
where rd^e = e^/(47reo/cBT'e) is the distance of closest 
approach for electrons. The electron-neutral collisional 
cross-section is Egn = ^ne = 1 x 10~^^ m^. The model 
used in Leake et al. 2012^^ neglected all but the first 
term on the right hand side of the generalized Ohm's law. 
Equation (pTj) . and used a parameter for the resistivity 
7^, rather than calculating it using Equation (p!2]) . 
Summary: 

The changes made to the model presented in Leake et al. 
2012^^, which have been implemented to make the physi- 
cal model more realistic, have a number of consequences. 
In the model in this paper we include charge exchange 
interactions, which, for the chromospheric plasma param- 
eters, effectively doubles the collisional coupling between 
ions and neutrals. Also, the ion-neutral scattering cross- 
section is larger here than in Leake et al. 2012^^, as we 
use a more realistic value, based on the theoretical work 
of Ref. Il6|. These two changes increase the neutral-ion 
collisional coupling, which reduces the neutral-ion mean 
free path by approximately an order of magnitude. As we 
will show in the following section, the Sweet-Parker width 
in this simulation is also smaller. As a result the ratio of 
Sweet-Parker width to neutral-ion collisional mean free 
path in the simulation in this paper is similar to that in 
Leake et al. 2012^^. 

In addition to the above changes, in the model used 
in this paper we evolve the system with a generalized 
Ohm's law where the resistivity is not a parameter of the 
simulations, but appropriately depends on the electron 
collision frequencies (uei^^en) and hence depends on the 
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local plasma parameters. This is more realistic than per- 
forming a parameter study over resistivity as was done 
in Leake et al. 20 12^^, and helps to confirm that the 
results found in that paper are valid for realistic chromo- 
spheric parameters. In addition to these differences the 
simulations in this paper have a magnetic Prandtl num- 
ber lower than the simulations of Leake et al. 2012^^, as 
we are using a different resistivity and more realistic ion 
and neutral transport coefficients. This will be discussed 
in more detail in the next section. 



B. Normalization 

The equations are non-dimensionalized by dividing 
each variable (C) by its normalizing value (Co). The set 
of equations requires a choice of three normalizing val- 
ues. Normalizing values for the length (Lq = 1 x 10^ m), 
number density (no = 3 x 10^^ m~^), and magnetic field 
(^0 = 2 X 10~^ T) are chosen. From these values the 
normalizing values for the velocity {vq = Bq/ ^yjaomino = 
2.52 X 10^ m.s"^), time (to = Lq/vq = 0.04 s), tempera- 
ture (To = B^/kB/J^ono = 7.69 x 10^ K), pressure (Po = 
B^/fiQ = 3.18 Pa), current density (jo = Po/(mo^o) = 
0.6 A.m^), and resistivity (?7o = /io^o'^o = 3.17 x 10^ Qm) 
can be derived. 



C. Initial conditions 

The simulation domain extends from -16Lo to 16Lo 
in the x direction and -4Lo to +4Lo in the y direction, 
with a periodic boundary condition in the x-direction 
and perfectly-conducting boundary conditions in the y- 
direction. The horizontal extent is chosen to be larger 
than the maximum horizontal extent of the current sheet 
length during the simulation (which is approximately 
6Lo). 

The initial neutral fluid number density is 200no (6 x 
2q18 ^-jr^^ ^Qj^ fluid number density is 0.2no 

(0.6 X 10^^ m~^). This gives a total (ion-|- neutral) mass 
density of 1.0 x 10~^ kg.m~^, and an initial ionization 
level {i/ji = ni/{ni + n^)) of 0.1 %. The initial electron, 
ion and neutral temperatures are set to 1.1 x 10~^To 
(8.46 X 10^ K). These initial conditions are consistent 
with lower to middle chromospheric conditions, based on 
ID semi-empirical models of the quiet Sun^^. For these 
plasma parameters, the isotropic neutral heat conduction 
is much faster than any of the anisotropic heat conduc- 
tion tensor components for either ions or electrons. This 
fact, coupled with the fact that the ion- neutral collision 
time is small compared with a typical dynamical time, 
means that thermal diffusion of the plasma is dominated 
by neutral heat conduction and is primarily isotropic. 

As in Leake et al. 2012^^, the initial magnetic config- 
uration is a Harris current sheet, 



Bo 



with an initial width of = 0.5Lo. Both the ionized 
pressure Pp = Pi -\- Pg and the neutral pressure P^ are 
increased in the current sheet in order to balance the 
Lorentz force of the magnetic field in the current sheet. 



Ppiy) = Pp 

Pn{y) = Pn 



1 F 

2 cosh^ {y/^^y 

1 Po-F 

2 cosh^ {y/^i>) ' 



(15) 
(16) 



where P = 2 x 10~^Po is chosen to maintain an approxi- 
mate ionization balance of 0.1% inside the current sheet. 

A relative velocity between ions and neutrals is there- 
fore required so that the collisional terms (both scatter- 
ing collisions and charge exchange) can couple the ion- 
ized and neutral fluids and keep the fluids individually in 
approximate force balance. 

For simplicity, using the fact that the ion-neutral drag 
due to charge exchange and scattering collisions are ap- 
proximately equal, we assume the following form for the 
initial ion velocity 



(P-1) tanh(7//A^) 
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A^cosh (y/X^) 
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(17) 



While this is not an exact steady state, the velocities that 
arise from the unbalanced forces in the initial condition 
are very small compared with the flows created by the 
onset of the reconnection process. 

To initiate the onset of the tearing mode instability, 
an additional, localized, small amplitude rotational flow 
perturbation is applied to the neutral and ion fluids, with 
the form 
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(19) 



Av;lncosh(?//A^), B = V A A^e^ (14) 



where vs^ = 5 x 10 "^vq and r = 

For these initial conditions, the background neutral- 
ion collisional mean free path, A^i = VT,n/^ni is 90 
m. The initial current sheet width is 0.5Po= 5 km, 
and so is much larger than A^i- Therefore we expect 
that the initial reconnection will behave like the classical 
single-fluid (i.e., when the neutrals and ions are coupled) 
Sweet-Parker model, and the initial current sheet will 
thin until it reaches the Sweet-Parker width, Ssp- To 
estimate t his Sweet -Parker width, we use the equation 
Ssp ^ ^/r]L/ /hqVa, use the initial resistivity rj, back- 
ground field strength and total number density (to get 
an estimate for va), and estimate that the reconnection 
region length L is equal to the spatial extent of the tear- 
ing perturbation {L = 5A^ = 2.5 x 10^ m). This gives 
an estimate of ~ 120 m, which is comparable to the 
background neutral-ion collisional mean free path, 90 m. 
If the current sheet thins until its width is comparable to 
the dsp^ its width will also be comparable to A^i, and so 
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the neutrals will begin to decouple from the ions. The de- 
coupling will result in the neutral pressure being unable 
to balance magnetic pressure and the sheet will thin fur- 
ther. Thus the current sheet width may decrease beyond 
the Sweet-Parker width, and eventually may become less 
then the mean free path. This scenario occurred in the 
simulations of Leake et al. 2012^^, although the physical 
model and initial parameters were different from the sim- 
ulation in this paper. In Leake et al. 2012^^, the back- 
ground neutral- ion collisional mean free path was 1350 
m, and the Sweet-Parker width had a range of 1330 m to 
8410 m (as the resistivity r] was varied as a parameter of 
the simulations). The current sheets in those simulations 
thinned to less than 1000 m, i.e., narrower than the mean 
free path. 



III. RESULTS 

Figure [1] shows the early evolution of the reconnection 
region resulting from the initial perturbation, on a sub- 
set of the simulation domain. The color indicates current 
density, and the solid lines show isocontours of A^, i.e., 
magnetic fieldlines. The y axis is stretched by a factor of 
40 to show the thin structures formed. The tearing in- 
stability grows, magnetic field reconnects at the X-point, 
and fiux is ejected. A Sweet-Parker-like reconnection re- 
gion is formed, where ions are pulled in by the magnetic 
field and drag the neutrals with them. 

The current sheet width and length as functions of 
time are shown in Figure [2l The width, Sgim^ is defined 
to be the half- width at half-max of the current sheet. 
The length, Lsim^ is defined as the x location along the 
line 7/ = at which the horizontal velocity of the ions 
reaches a maximum (i.e., the location of maximum out- 
fiow). These two definitions are the same as was used in 
Leake et al. 2012^^. As can be seen in Figure [21 the cur- 
rent sheet width falls from a value of 0.037Lo = 370 m at 
time t = 800^0 to a value of 3.3 x 10~^Lo = 33 m by the 
end of the simulation. This final width is less than the 
neutral- ion collisional mean free path at this time (95 m). 
Thus during the simulation the current sheet has under- 
gone a transition from being thicker than the mean free 
path to being thinner than it, as predicted. The current 
sheet also lengthens after t = 1500to, increasing to 2.9Lo. 
The aspect ratio {Ssim/ Lsim) at t = 1800to is 1.1 x 10~^. 

Figure [3] shows profiles across the current sheet at 
X = of the ion density, current density, vertical ion fiow, 
and the difference between the vertical neutral fiow and 
the vertical ion fiow, at three different times in the simu- 
lation (1408.9^0, 1661.9^0, and 1803.5to). As the tearing 
instability develops. Panels (a) and (b) show the current 
density and ion density building into a sharp structure 
with the same gradient scale. Note that at = 0, the 
ionized fraction {ni/{ni + n^)) rises by an order of mag- 
nitude from 0.1% at t = to 1.2% at the end of the sim- 
ulation {t = 1807to), though the plasma remains weakly 
ionized throughout the domain. Panels (c) and (d) show 



an increase in the ion infiow and an increase in the differ- 
ence between neutral and ion infiows as the current layer 
develops. The amplitude of the ion infiow reaches a max- 
imum of 1.5 X 10~^^'o5 while the neutral infiow is much 
slower than the ion infiow. The peak difference between 
neutral and ion infiows at t = 1803.5 s is 1.4 x 10~^vo, 
which means that the peak neutral infiow is 1 x 10 "^^o, 
less than a tenth of the peak ion infiow. Clearly the neu- 
trals have become decoupled from the ions, and are being 
pulled into the reconnection region much slower than the 
ions. 

Figure S] shows profiles of the ion density, current den- 
sity, ion outflow and the difference between the neutral 
and ion outflows, within the current sheet, along the line 
y = at three different times. The wave-like oscillations 
seen in the current density proflle between x = and 
X = I/O in Panel (b) are caused by the onset of the sec- 
ondary tearing instability, which will be discussed later 
in this section. Also note the increase in ion density and 
moderate increase in current density magnitude outside 
the reconnection region, x > 3.5Lo, at later times. This 
behavior is characteristic of the ID magnetic fleld recon- 
nection in a weakly ionized plasma considered by Ref. ITsI 
and Ref. |l9|. As discussed in the above referenced stud- 
ies, such a recombination facilitated process can lead to 
resistivity independent reconnection for a limited range 
of plasma parameters. However, as will be shown here 
and as was shown by Leake et al. 2012i^, by allowing 
for plasma outflow in a 2D conflguration and with given 
plasma conditions, the reconnection rate can be acceler- 
ated beyond the ID result. 

Figure m Panels, (c) and (d), show that the neutral 
and ion outflow are very well coupled with the difference 
between neutral and ion outflow being approximately 4 
orders of magnitude smaller than the ion outflow. To 
highlight the fact that the neutral inflow is decoupled 
from the ion inflow, but that the outflows are coupled. 
Figure [5fa) shows a subset of the simulation domain at 
t = 1661.9^0- The solid lines show isocontours of A^. The 
top left quadrant shows the ion inflow velocity, and the 
top right quadrant shows the neutral inflow velocity. At 
this time, the current sheet width 6sim = 0.0045Lo = 45 
m, which is smaller than the neutral-ion collisional mean 
free path A^i = 0.0095Lo = 95 m. The ions have a larger 
maximum amplitude of inflow than the neutrals, which is 
a direct consequence of the fact that the current sheet is 
now narrower than the neutral-ion collisional mean free 
path and the neutrals are not dragged into the current 
sheet as quickly as the ions. Figure [5fb) shows the out- 
flow velocities. The outflows are well coupled in the sense 
that the difference between ion and neutral outflow is 
negligible compared to the magnitude of the ion outflow. 
This situation, where the inflow is decoupled but the out- 
flow is coupled, was also observed in the earlier simula- 
tions by Leake et al. 2012^^, and is a consequence of the 
neutral-ion mean free path lying between the inflow and 
outflow scales of the current sheet, Ssim < ^ni < Lgim- 

The maximum outflow speed during the simulation is 
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0.0126'Uo = 3 km/s which, although smah compared to 
typicahy measured spicule flows (up to 150 km/s), is close 
to the flows observed in chromospheric jets (~ 10 km/s). 
The Alfven speed just upstream of the reconnection cur- 
rent sheet is 0.006'Uo, and so the outflow can be up to 
twice the local upstream Alfven speed. As the Alfven 
speed scales with magnetic fleld strength and density, it 
is possible that larger outflows can be obtained by using 
either larger flelds {Bq = 20 G here), corresponding to 
more active regions of the Sun, or by using lower num- 
ber densities, corresponding to conditions in the upper 
chromosphere. 

Figure [6] displays contributing terms to the y com- 
ponent of the total momentum equation at time t = 
1661.9^0- The top left quadrant shows the vertical gra- 
dient of the ionized (ion+electron) pressure, and the top 
right quadrant shows the vertical gradient in the neutral 
pressure. The bottom half shows the y component of the 
Lorentz force. The color scale highlights that on the scale 
of Ani, it is the ionized pressure that balances the Lorentz 
force (in the y direction), and on scales larger than A^i 
it is the neutral pressure that balance the Lorentz force. 

The separation of inflow velocities has consequences for 
the reconnection occurring in the current sheet. If ions 
are pulled in faster than neutrals, with equal outflows, 
then there will be an excess of ions above the ionization 
equilibrium, and recombination will occur. In Leake et 
al. 20 li*^^ we found that recombination was as efficient at 
removing ions from the reconnection region as the out- 
flow caused by reconnecting fleld. Figure [7] shows the 
four components contributing to in the ion continu- 
ity equation for the simulation described here. The top 
left quadrant shows the loss due to recombination, the 
bottom left quadrant shows the loss due to outflow (hor- 
izontal gradient in horizontal momentum of ions), the 
top right quadrant shows the gain due to inflow (vertical 
gradient in vertical momentum), and the bottom right 
quadrant shows the gain due to ionization. Unlike in 
Leake et al. 2012^^, here we observe that, within the re- 
connection region, outflow is larger than recombination, 
with a negligible contribution from ionization. An im- 
portant distinction between the simulation in this paper 
and those in Leake et al. 2012^^, is that in the present 
simulation the current sheet is never in steady state. This 
is indicated by the monotonic decrease in current sheet 
width throughout the simulation, as shown in Figure [21 
Panel (a), as well as the temporal evolution of the cur- 
rent sheet, shown in Figure [51 Hence the steady state 
analysis, which leads to the reconnection rate being de- 
termined by a balance between outflow, ionization and 
recombination, is not strictly valid. 

We now look at the temporal evolution of the recon- 
nection rate in the current sheet. The normalized recon- 
nection rate is deflned by 



V Jr. 



(20) 



current density, j^, located at {x^y) = (xj,0). Bup is 
Bx evaluated at {xj^5sim)^ where Ssim has already been 
deflned as the half- width half-max in out of plane current 
density, is the relevant Alfven velocity deflned using 
Bup and the total number density (n*) at the location of 
jmax- To obtain a local estimate of resistivity we use 



J max 



(21) 



with Rg* +Rg"^ evaluated at the location of jmax- The six 



quantities r^*, j„ 



Bup, v% and M,,ir, 



are shown as 



functions of time in Figure [H The resistivity 77* has two 
contributions, one from electron- ion collisions and the 
other from electron- neutral collisions. Taking the ratio 
of these two terms gives 



(22) 



where Sgn = 1-16 x 10~^^ m^, and depends on elec- 
tron temperature as was described in Section 2.1. The 
ion, electron, and neutral temperatures have been as- 
sumed equal to each other in this estimate. For the ini- 
tial temperature of 8460 K, Sgi = 1.2 x 10~^^ m^, and at 
t = 0, Un/rii = 10^. This gives a ratio of ^ = 1.5. 
At t = 1780^0, the electron temperature is 9230 K, 
Sei = 1 X 10~^^ m^ and n^/n^ = 80, which gives a ratio of 
^ = 0.13. Thus initially the electron- neutral contribu- 
tion is of the same order as the electron-ion contribution, 
but toward the end of the simulation, the electron-ion 
contribution dominates due to the increase in ion den- 
sity. 

The increase of the current density with time is primar- 
ily responsible for the monotonic increase of Mgimi from 
0.001 at t = 400^0 to 0.056 at the end of the simulation. 
The Sweet-Parker model for magnetic reconnection pre- 
dicts the normalized reconnection rate M = 1/ \/~S where 
S is the Lundquist number of the reconnecting current 
sheet. The Lundquist number in this simulation, Ssirm 
is deflned by 



Ss 



7)* /■* 



(23) 



At t = 1750 s, before the plasmoid instability has de- 
veloped, 77* = 1.205 X 10~^7/o, ^sim = 2.5Lo, and 



0.0066^0, which gives a value of Sg^ 



1.37 X 10"^. 



Here, jmax is the maximum value of the out of plane 



Hence the Swee t-Parker m odel would predict a reconnec- 
tion rate of 1/ Vl.7 x 10^ = 0.0085. In this simulation we 
see the reconnection rate rise from below 0.001 to above 
0.048 at t = 1750 s. As the current sheet width be- 
comes less than A^i and the neutrals decouple from the 
ions, it is clear that a single-fluid Sweet-Parker theory 
is no longer applicable. The reconnection rates observed 
in this simulation are sufficient to explain the observed 
lifetimes and reconnection rates of spicules, but not chro- 
mospheric jets. 

At t = 1790^0, the current sheet becomes unstable to 
the secondary tearing mode, known as the plasmoid in- 
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stability (e.g., Ref. '20'). Figure [9] shows the ion and cur- 
rent density, and selected isocontours of A^. A plasmoid 
can be seen at the origin at t = 1803.5to. The laminar 
current sheet can also be seen to break up into thinner 
current sheets. In the context of a single-fluid theory, Ref. 
[20! predicted that the onset of the plasmoid instability oc- 
curs when the aspect ratio decreases to 1/200 and this 
was the case in the multi-fluid simulations of Leake et al. 
2012i^. At t = 1700^0 the aspect ratio asim = ^sim/Lsim 
already is 1/500, yet the plasmoids are not observed until 
t = 1790^0 when the aspect ratio has fallen further still. 
At this time we do not have a thorough understanding 
of why the onset of the plasmoid instability is different 
between the present simulation and our previous work. 
However, we note that the magnetic Prandtl number, de- 
fined by the ratio of viscous to resistive diffusivity, 

Pr = M^, (24) 

in this simulation is approximately 0.4. This is at least 
an order of magnitude lower than that used in Leake et 
al. 2012^^, where Pr ranged from 5 for high resistivity to 
200 for low resistivity. The lower Pr in the simulation in 
this paper means that viscous diffusivity is less important 
and this allows for stronger shearing of the reconnection 
outflows that may lead to stabilization of the secondary 
tearing instability^^. A detailed investigation of the sta- 
bility of weakly ionized reconnection current sheets to 
secondary tearing instability is left to future work. 

Fi gure 1101 shows the reconnection rate MgiYn and jmax 
in the interval 

t = [1650, 1803. 5]to, i.e. a zoom in of Figure [51 Panels 
(b) and (f). At t = 1790to, there is a change in the time 
derivative of jmax as the plasmoids create thinner current 
sheets and larger currents. This creates an increase in the 
reconnection rate. The further study of the non-linear 
evolution of the plasmoids is limited in this investigation, 
as the current sheet width approaches the limit of the 
resolution of the numerical simulation. 

In this paper magnetic reconnection in the solar chro- 
mosphere is simulated using a partially ionized reacting 
multi-fluid plasma model, which is able to capture the 
physical effects of decoupling between the ion and neu- 
tral fluids, occurring when scales associated with the re- 
connection region become comparable to the neutral-ion 
collisional mean free path. The number densities, mag- 
netic field strengths, ionization levels and temperature 
are consistent with lower to middle chromospheric con- 
ditions in a quiet Sun region. 

These simulations improve upon the earlier work of 
Leake et al. 2012^^ by including the effects of charge- 
exchange interactions, and using more realistic values for 
the cross-sections of scattering collisions and transport 
coefficients. In addition, the resistivity r] depends on the 
local plasma parameters, rather than being a fixed value. 
Even so, our simulation lies in the same regime as those 
of Leake et al. 2012^^. This regime consists of an initial 
current sheet width Ssim which is much larger than the 



Sweet-Parker width 5sp^ which itself is comparable to 
the neutral-ion collisional mean free path: 

Ssim{t=0) > ^SP ^ Ani- (25) 

In this regime, the current sheet width tends towards 
the Sweet-Parker width Ssp (and hence towards Xni)- 
When the current sheet width becomes comparable to 
the neutral-ion mean free path, decoupling of neutrals 
from ions allows the sheet width to fall below both the 
Sweet-Parker width and the mean free path: 

(26) 

Once the current sheet width is smaller than the neutral- 
ion collisional mean free path, the steady state recon- 
nection rate is determined by a combination of ion out- 
flow and recombination. In the simulations of Leake et 
al. 2012i^, recombination was found to be as important 
as outflow at removing ions from the reconnection re- 
gion, and steady state reconnection rates faster than the 
Sweet-Parker prediction were seen. On the other hand, 
in the simulation presented in this paper no steady state 
is achieved. Instead we find an instantaneous reconnec- 
tion rate which is faster than the Sweet-Parker prediction 
due primarily to the separation of the current sheet width 
and the neutral-ion collisional scale and the resulting or- 
der of magnitude enhancement of the ion density in the 
current sheet with respect to the upstream value. 

An important characteristic of the solar chromosphere 
is that it has a high degree of spatial variation of the 
plasma temperature and density due to stratification and 
horizontal structure, and also a high degree of spatial 
variation of the small scale magnetic field. This situation 
does not easily lend itself to investigating phenomena in a 
"typical" chromosphere. In particular, it is not clear that 
the regime we find in this simulation, where 5sim < ^ni^ 
is a ubiquitous occurrence in the chromosphere. To ad- 
dress this issue, we can estimate the altitude profiles of 
both the Sweet-Parker width and neutral-ion collisional 
mean free path using a semi-empirical ID model of the 
quiet Sun^^, and a few assumptions. The neutral-ion col- 
lisional mean free path depends on temperature and ion 
density via the thermal speed and neutral-ion collisional 
frequency. The temperature and density are provided 
by the ID atmospheric model of Ref. [22|. To obtain an 
altitude profile for the Sweet-Parker width we can re- 
place the characteristic system size L with Ssp/crsp in 
the equation Ssp ^ ^yrjL/vAl^o to obtain 

Ssp ^ ^ . (27) 

crspVA/i^o 

Taking an aspect ratio at which the plasmoid instabil- 
ity should theoretically set in gives a lower bound on the 
Sweet-Parker aspect ratio asp and hence an upper bound 
on the Sweet-Parker width Ssp- This critical value varies 
in the literature, and so we take two extremes, 1/1000 
(of the order of the smallest value seen in the simula- 
tion in this paper), and 1/50 (values observed in high 
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plasma-beta simulations by Ref. ^2§). An estimate of 
the magnetic field strength in the chromosphere is also 
required to estimate 6s p in the chromosphere. We use 
the range [5,1000] G. Using these two sets of extremes 
(one for asp and one for magnetic field strength) gives a 
range of altitude profiles for Ss p - Figure [TT] depicts the 
two extreme profiles of ^^p, as well as the altitude profile 
of the neutral-ion collisional mean free path. This plot 
shows that for the ranges of aspect ratios and magnetic 
field strengths considered here, the Sweet-Parker width 
can be greater than, but comparable to, the neutral-ion 
collisional mean free path in a height range of 650 km 
to 2000 km above the solar surface. This essentially en- 
compasses the entire chromosphere. Hence, the regime 
that our simulation self-consistently produces, where the 
reconnection transitions from being fully coupled to de- 
coupled, is possible for a wide range of chromospheric 
conditions. This transition from coupled to decoupled 
reconnection is accompanied by an increase in the recon- 
nection rate. 

Recombination-aided fast reconnection in a weakly 
ionized plasma that is out of ionization balance was con- 
sidered by Ref. [H and Ref. 19. Those theoretical in- 
vestigations modeled the reconnection region as ID and 
time independent with a width which was already less 
than the neutral-ion collisional mean free path (denoted 
Lad in Ref. [lH, and denoted A^i in this paper). Hence 
the reconnection region was assumed to already have an 
excess of ions due to the decoupling of neutrals and ions. 
The simulation performed in this paper shows that for 
initial conditions relevant to the solar chromosphere, and 
a current sheet initially thicker than A^i, the system can 
self-consistently evolve to a current sheet width thinner 
than Xni. However, in this simulation, the recombina- 
tion is less important than the outfiow of ions, whereas 
in the ID models of Ref. 18 and Ref. 0, the recombina- 
tion was assumed to be the only sink for the ions in the 
reconnection region. 

As in Leake et al. 2012^^, we measure "fast" recon- 
nection rates close to M = 0.1. Generally, fast reconnec- 
tion requires either Hall, kinetic, or localized resistivity 
effects^^, or secondary tearing^^^^^. Hall effects are neg- 
ligible in this simulation as the ion inertial scale is less 
than the smallest current sheet width in this simulation. 
In addition, the onset of the secondary tearing mode only 
increases the reconnection rate after t = 1750 by approx- 
imately 15%. In this simulation we obtain normalized 
reconnection rates above 0.05 without either Hall effects 
or secondary tearing. Given that we also have neither a 
localized increase in resistivity in the reconnection site, 
nor kinetic effects, we can conclude that fast reconnec- 
tion is obtained solely due to the decoupling of neutrals 
from ions when the current sheet width is less than the 
neutral-ion collisional mean free path. Note, however, 
that Figure [TT] shows that the Sweet-Parker width can 
be as low as 1 cm near the top of the chromosphere, and 
hence potentially smaller than the ion inertial scale, and 
so it is possible for Hall effects to be important in some 



regions of the chromosphere. 

The values of reconnection rate in this simulation 
are comparable to estimated rates from observations of 
spicules. This fact along with the fact that we see out- 
fiows of 3 km/s, comparable to observations of chromo- 
spheric jets, provides strong evidence that weakly ionized 
magnetic reconnection is capable of explaining the ubiq- 
uitous observance of transient phenomena in the chromo- 
sphere. The plasmoids observed in this simulation, when 
the current sheet becomes unstable to the secondary tear- 
ing mode (plasmoid instability), may also explain the ob- 
servations of "blobs" of plasma within chromospheric jet 
outfiows^^^ This paper shows that non-equilibrium two 
fiuid (ion + neutral) effects can be important for recon- 
nection on scales larger than those at which Hall and 
kinetic effects become important in the chromosphere. 
Hence it is vital that studies of chromospheric reconnec- 
tion, and the subsequent dynamics and heating this re- 
connection causes, include ion-neutral collisional physics, 
as has been done in the model presented in this paper. 
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FIG. 1. Formation of reconnect ion region. Current den- 
sity (jz/jo) is shown by color scale, and 5 contour values of 
A^, [-0.02, -0.0175, -0.015, -0.0125, -0.01]BoLo, show rep- 
resentative fieldlines. Note that the ^/-coordinate is expanded 
by a factor of 40. 
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FIG. 2. Width and length of the current sheet as a function of 
time. The width is defined by the half-width half- maximum 
based on the current density jz . The length is defined as the 
distance from x = to the x location of maximum ion outflow 
Vi^y on the y — line. The current sheet width monotonically 
decreases during the simulation after t=800 s, and eventu- 
ally becomes less than the background neutral-ion collisional 
mean free path Ani, which is 0.0095Lo = 95 m at the end of 
the simulation. Before 800 s, the tearing instability has not 
developed enough to give a uniquely defined current sheet 
width and length. As the non-linear reconnection continues, 
the current sheet is seen to elongate after t=1500 s. 
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FIG. 3. Profiles across the current sheet (x = 0). Panel 
(a) shows the ion density (rii/no), Panel (b) shows the cur- 
rent density (jz/jo), Panel (c) shows the ion inflow (vi^y/vo), 
and Panel (d) shows the difference in neutral and ion inflow 
,y — Vn,y)/vo). This figurc demonstrates the sharpening 
in the current density and ion density as the current sheet 
collapses. It also indicates that the peak in \vn,y \ is approxi- 
mately 0.1 times the peak in \vi^y\, i.e., the inflow is decoupled. 
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FIG. 4. Profiles along the current sheet = 0). Panel (a) 
shows the ion density (n^/no). Panel (b) shows the current 
density (jz/jo)- Panel (c) shows the ion outflow (vi^x/vo). 
Panel (d) shows the difference in neutral and ion outflow 
{{vi,x — Vn,x)/vo). The difference in ion and neutral outflow 
is small compared to the ion outflow, and so the neutrals and 
ions fluids are coupled. The x location of the peak in the ion 
outflow can be seen to move outward in time (Panel (c)). 
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FIG. 5. Illustration of the decoupling of ion and neutral 
inflow during reconnection. At t = 1661. 9to, the current 
sheet width 5sim is less than the background neutral-ion col- 
lisional mean free path \ni. Panel (a): The ion (top left 
quadrant) and neutral (top right quadrant) inflows {vi,y/vo 
and Vn,y/vo). The solid lines show 5 contour values of Az, 
[-0.01, -0.00975, -0.0095, -0.00925, -0.009]BoLo. The vec- 
tors show ion flow on the bottom left quadrant and neutral 
flow on the bottom right quadrant. Panel (b): The ion and 
neutral outflows {vi^x/vo and Vn,x/vo): The contour lines and 
vectors are the same as in the top panel. 
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FIG. 6. Contributions to the plasma and neutral mo- 

dP 

mentum equations. The top left quadrant shows where 
Pp — {Pi + Pe), and the top right quadrant shows 
The bottom half shows {j A B)y. The color is on a log- 
scale. This shows that balances {j A B)^ on the col- 

lisional scale Xni, while balances {jAB)^ on scales 

greater than Xni. The solid lines are 5 contour values of Az, 
[-0.01, -0.00975, -0.0095, -0.00925, -0m9]BoLo. 
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FIG. 7. The steady state reconnection region showing con- 
tributing sources and sinks of ions in the current sheet 
(in units of L^^t^^): Top left quadrant shows rate of 
loss of ions due to recombination. Bottom left quadrant 
shows rate of loss of ions due to outflow Top 
right quadrant shows rate of gain of ions due to inflow 

Bottom right quadrant shows rate of gain of ions 

due to ionization. The solid lines are 5 contour values of 
A^, [-0.01,-0.00975, -0.0095, -0.00925, -0.009]5oLo. This 
shows that ionization is negligible, and that outflow is larger 
than recombination. 



400 800 1200 1600 



400 800 1200 1600 




0.00 



400 



800 



1200 1600 



310 



305 



g 300 



+ 



295 



290 



285 




400 



800 



1200 1600 



0.002 
0.000 



T 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1— 




_L_L_ 



400 800 1200 1600 




400 



800 



1200 1600 

*/*0 



FIG. 8. Temporal evolution of the reconnect ion rate. Panel 
(a): Effective resistivity ?7*/?7o. Panel (b): jmax/jo, the max- 
imum value of the current density, located at {x^y) — (xj, 0), 
within the reconnection region. Panel (c): Bup/Bo where Bup 
is evaluated at (xj^Ssim) and Ssim is the half- width at half- 
max of the current sheet. Panel (d): n*/no, the total density 
inside the current sheet at (xj,0). Panel (e): v\/vo, the cur- 
rent sheet Alfven speed. Panel (f): Msim — r/^jmax/v^Bup. 
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FIG. 9. Formation of plasmoids. Current density (jz/jo) on 
the right, and ion density rii/no on the left. The black line 
is the Az — — 0.01245^0 Lo contour and the white line is the 
Az = —0. 01237^0 Lo contour. Note that the y-coordinate is 
expanded by a factor of 40. 
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FIG. 10. Temporal evolution of the reconnection rate dur- 
ing the latter stages when the plasmoid instability sets in. 
Left panel: jmax/jo, the maximum value of the out of plane 
current density. Right panel: Mgim — T]jmax/v\Bup, the nor- 
malized reconnection rate. 
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FIG. 11. Estimates of the Sweet-Parker width {6sp) and neu- 
tral ion colhsional mean free path (Xni) using the FALC model 
of the solar chromosphere^^. The solid lines show two ex- 
tremes of the calculations of the resistive length. The higher 
value uses an aspect ratio of a = 1/1000 and a field strength 
of 5 G, and the smaller uses a current sheet aspect ratio of 
cr = 1/50 and a field strength of 1000 G. The dot-dashed line 
is the neutral- ion collisional mean free path, Xni. This figure 
shows that for a range of heights in the chromosphere, the 
Sweet-Parker width can be larger than, but comparable to, 
the neutral-ion collisional mean free path, Xni. 



